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The initial stages of planet formation in circumstellar gas discs proceed via dust grains that 
collide and build up larger and larger bodies 1 . How this process continues from metre-sized 
boulders to kilometre-scale planetesimals is a major unsolved problem 2 : boulders stick to- 
gether poorly 3 , and spiral into the protostar in a few hundred orbits due to a head wind from 
the slower rotating gas 4 . Gravitational collapse of the solid component has been suggested 
to overcome this barrier 1 ' 5 ' 6 . Even low levels of turbulence, however, inhibit sedimentation 
of solids to a sufficiently dense midplane layer 2 7 , but turbulence must be present to explain 
observed gas accretion in protostellar discs 8 . Here we report the discovery of efficient grav- 
itational collapse of boulders in locally overdense regions in the midplane. The boulders 
concentrate initially in transient high pressures in the turbulent gas 9 , and these concentra- 
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tions are augmented a further order of magnitude by a streaming instability 10-12 driven by 
the relative flow of gas and solids. We find that gravitationally bound clusters form with 
masses comparable to dwarf planets and containing a distribution of boulder sizes. Gravita- 
tional collapse happens much faster than radial drift, offering a possible path to planetesimal 
formation in accreting circumstellar discs. 

Planet formation models typically treat turbulence as a diffusive process that opposes the 
gravitational sedimentation of solids to a high density midplane layer in circumstellar discs 7 - 13 . Re- 
cent models of solids moving in turbulent gas reveal that the turbulent motions not only mix them, 
but also concentrate metre- sized boulders in the transient gas overdensities 9 formed in magnetoro- 
tational turbulence 14 , in giant gaseous vortices 15 ' 16 , and in spiral arms of self-gravitating discs 17 . 
Short-lived eddies at the dissipation scale of forced turbulence concentrate smaller millimetre-sized 
solids 18 . 

Some simulations mentioned above 9 ' n ' 12 were performed with the Pencil Code, which solves 
the magnetohydrodynamic (MHD) equations on a three-dimensional grid for a gas that interacts 
through drag forces with boulders. Boulders are represented as superparticles with independent 
positions and velocities, each having the mass of a huge number of boulders but the aerodynamic 
behaviour of a single boulder. We have now further developed the Pencil Code to include a fully 
parallel solver for the gravitational potential of the particles (see Supplementary Information). The 
particle density is mapped on the grid using the Triangular Shaped Cloud assignment scheme 19 
and the gravitational potential of the solids is found using a Fast Fourier Transform method 20 . 
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This allows us, for the first time, to simulate the dynamics of self-gravitating solid particles in 
magnetised, three-dimensional turbulence. 

We model a corotating, local box with linearised Keplerian shear that straddles the protoplan- 
etary disc midplane and orbits the young star at a fixed distance. Periodic boundary conditions are 
applied. An isothermal equation of state is used for the gas, while the induction equation is solved 
under the ideal MHD assumption of high conductivity. Magnetorotational instability 14 drives tur- 
bulence in Keplerian discs with sufficient ionisation 21 , producing in our unstratified models turbu- 
lence with Mach number Ma ks 0.05 and viscosity a xs 10~ 3 , a realistic value to explain observed 
accretion rates 8 . The ionisation fraction in the dense midplanes of protoplanetary discs may be 
insufficient for the gas to couple with the magnetic field to drive magnetorotational instability 21 . 
In the Supplementary Information we therefore describe unmagnetised models as well. 

Solid objects orbit the protostar with Keplerian velocity v K in the absence of gas drag. A 
radial pressure gradient partly supports the gas, however, so it orbits at sub-Keplerian velocity, 
with Av = v g — vk < 0. As a result, large (approximately metre-sized) solid objects feel a strong 
head wind that causes them to drift radially inwards 4 with a maximum drift velocity Av. They 
also feel gas drag as they fall toward the disc midplane in the effective gravity field of the star. 
A sedimentary midplane layer forms with a width determined by a balance between settling and 
turbulent diffusion 7 ' 13 . 

We present three types of models: (1) without self-gravity, with 128 3 zones and 2 x 10 6 
particles, run for 100 orbits, to study the interplay between the streaming instability 10 and concen- 
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tration by transient high pressures; (2) with self-gravity and boulder collisions, with 256 3 zones and 
8 x 10 6 particles run for 27 orbits, to study gravitational collapse; and (3) models with self-gravity 
but no magnetorotational turbulence (presented in the Supplementary Information). Magnetorota- 
tional turbulence is given 10 orbits to reach steady state before we turn on drag force and vertical 
gravity, to avoid any influence of the initial conditions on the sedimented midplane layer. We fix 
the global solids-to-gas bulk density ratio at the canonical galactic value of e = 0.01, but two 
values of the radial drift are considered: low drift with Av = — 0.02c s , where c s is the isothermal 
sound speed, and moderate drift with Av = — 0.05c s , depending on the assumed radial pressure 
support (values up to Av = —[0.2 . . . 0.5]c s are possible 13 , but are not considered here). 

For the simulations without self-gravity we consider a fixed particle size parameterised by the 
dimensionless friction time Q-^Ti = 1.0, where i?K is the local Keplerian rotation frequency and T{ 
is the time- scale over which gas and solids reach equal velocity. At an orbital distance r = 5 AU 
this corresponds to boulders of approximately one metre in diameter. Figure [T] shows the space- 
time topography of the sedimented midplane layer. The streaming instability increases the density 
of boulders in regions where they have already been concentrated by transient high pressures 9 . 
Increasing radial pressure support from Av = — 0.02c s to — 0.05c s reduces the concentration by 
streaming instability, although the local solids-to-gas density ratio still reaches 200. 

Gravitational collapse of discrete solid objects produces virialised clusters unable to contract 
further 22 in the absence of mechanisms to dynamically cool the cluster-that is, reduce the local 
rms speed. Two processes that we consider can be important: drag force cooling and collisional 
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cooling. Drag force cooling occurs because part of the kinetic energy exchanged between the par- 
ticles and the gas is dissipated. Collisional cooling is produced by the highly inelastic collisions 
between boulders, transferring kinetic energy to heat and deformation. Collisional cooling occurs 
generally in simulations of resolved collisions in planetary rings 23 . In the Supplementary Informa- 
tion we describe how we treat collisional cooling numerically in the self-gravitating simulations by 
damping the rms speed of the particles in each grid cell on a collisional time-scale. We have found 
that in the absence of collisional cooling, gravitational collapse still proceeds if the total surface 
density (of solids and gas) is augmented by 50%. Collisional cooling is thus not a prerequisite 
of the collapse, but does allow it to occur in somewhat less massive discs. We ignore all other 
effects of the collisions, such as coagulation and collisional fragmentation. Collisional cooling and 
self-gravity are turned on after 20 orbits in the self-gravitating simulations. 

Our chosen scale-height-to-radius ratio of H/r = 0.04 gives a gas temperature of T = 
80 K at an orbital radius of r = 5 AU. We choose for the 256 3 self-gravitating run the uniform 
gas volume density to be consistent with the midplane of a disc with surface density of E gas = 
300 g cm~ 2 . This corresponds to approximately twice the minimum mass solar nebula (MMSN) 
at 5 AU from the (proto-)Sun. An alternative theory for giant planet formation, the disc instability 
hypothesis 24 - 25 , requires column densities at least 20 times higher than the MMSN for gravitational 
fragmentation of the gaseous component of the disc to occur. 

We have examined the numerical convergence of our models with resolutions ranging from 
64 3 to 256 3 zones (see Supplementary Information). The peak particle density on the grid increases 
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with increasing resolution, because of less smoothing in the particle-mesh scheme at higher reso- 
lution, resulting in a decrease in the column density threshold for gravitational collapse. Although 
we have not yet fully converged, our results appear to provide good upper limits to the column 
density for which collapse can occur. For the self-gravitating simulation we consider boulders 
with friction times distributed among Qkt~{ = 0.25, 0.50, 0.75, 1.00. At r = 5 AU in our chosen 
disc model, these correspond to radii of 15-60 cm. Consideration of multiple boulder sizes is vital 
since differential aerodynamic behaviour could inhibit gravitational instabilities 26 . The size range 
covers roughly half of the two orders of magnitude in particle radius produced by coagulation of 
microscopic grains 27 . Smaller particles are ignored since they are unlikely to separate from the gas 
and participate in gravitational collapse. In case of widespread collisional fragmentation, e.g. in 
the warmer terrestrial planet formation region, up to 80% of the solid material may be bound in 
small fragments 27 ' 28 , in which case we must implicitly assume an augmentation in solids-to-gas 
ratio of up to 5. 

In our self-gravitating model we set Av = — 0.02c s , but show in the Supplementary Informa- 
tion that gravitationally bound clusters also form for Av = — 0.05c s , with a factor of two increase 
in column density threshold. The Supplementary Information also documents that typical boulder 
collisions happen at speeds below the expected destruction threshold 3 . We caution, however, that 
material properties, and thus destruction thresholds, of the boulders are poorly known. Higher 
resolution studies, and an improved analytical theory of collision speeds that takes into account 
epicyclic motion, will be needed to determine whether collision speeds have converged, given an 
unexplained factor 3 difference for = 1 particles between typical relative speeds within cells 
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(si 5 ms ) and the expected collision speed of well-mixed particles. 

The development of gravitational instability in the 256 3 run is shown in Figure [2] The four 
different boulder sizes have accumulated in the same regions {central panel) already before self- 
gravity is turned on, demonstrating that differential drift does not prevent density enhancement 
by streaming instability. Gravitationally bound boulder clusters form, with a Hill radius — within 
which the gravity of the cluster dominates over tidal forces from the central star — that increases 
steadily with time {inserts) as the clusters accrete boulders from the surrounding flow. 

We show in Figure|3]the peak density and the mass of the most massive gravitationally bound 
cluster as a function of time. The cluster consists of particles of all four sizes, demonstrating 
that different boulder sizes can indeed take part in the same gravitational collapse, despite their 
different aerodynamical properties and drift behaviour. At the end of the simulation the most 
massive cluster contains 3.5 times the mass of the dwarf planet Ceres. The cluster mass agrees 
roughly with standard estimates from linear gravitational instability 5 , when applied at r = 5 AU to 
the locally enhanced column densities (see Supplementary Information). 

Models lacking magnetic fields, and thus magnetorotational turbulence, are described in the 
Supplementary Information. Here sedimentation occurs unhindered until the onset of Kelvin- 
Helmholtz instabilities driven by the vertical shear of gas velocity above the midplane boulder 
layer. Strong boulder density enhancements in the mid-plane layer nevertheless still form, al- 
though for moderate drift an increase of the solids-to-gas ratio from 0.01 to 0.03 (perhaps possible 
due to radial variation in boulder drift speeds 6 or photoevaporation of the gas 29 ) was needed to 
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obtain strong clumping. Magnetorotational turbulence thus has a positive effect on the mid-plane 
layer's ability to gravitationally collapse, although collapse can occur without it as well. 

The Supplementary Information also includes a model with an adiabatic equation of state 
and explicit gas heating due to energy dissipated by drag and inelastic collisions. We find that gas 
heating does not prevent collapse. The maximum temperature reached is not even high enough to 
melt ice, although that may change with the formation of massive bodies with escape velocity near 
the sound speed. 

Our proposed path to planetesimal formation depends crucially on the existence of a dense 
sedimentary layer of boulders. Future investigations should focus on the formation and survival of 
such layers in light of processes like coagulation, collisional fragmentation and erosion 28 . Espe- 
cially important are higher resolution studies of collision speeds and an improved analytical theory 
of collisions that includes the epicyclic motion of particles. 
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Figure 1 Topography of the sedimented particle layer in models without self-gravity 
or collisional cooling, a) The azimuthally averaged vertical column density E v of metre- 
sized boulders (with f2 K r { = 1) as a function of radial coordinate x and time t, in a model 
where the particles feel gas drag, but the gas does not feel drag from the particles. Radial 
drift is evident from the tilted bands (particles crossing the inner boundary reappear at 
the outer). Transient regions of mildly increased gas pressure temporarily concentrate 
boulders. The gas orbits slightly slower on the outer edge of these high pressure regions 
and slightly faster on the inner edge, resulting in a differential head wind that forces boul- 
ders towards their centres 9 30 , b) Including the drag force from the particles on the gas 
allows for the development of the streaming instability, seeded by the existing radial den- 
sity enhancements. The streaming instability occurs where the collective drag force of 
the solids forces the gas to locally move with an orbital speed that is closer to Keplerian, 
reducing the gaseous head wind that otherwise causes boulders to drift radially. Solids 
then drift into already overdense regions from further out, causing runaway growth in the 
local bulk density of solids, c) The column density when the radial pressure support is 
increased from Av = -0.02c s to -0.05c s . Radial density enhancements become narrower 
and shorter-lived due to downstream erosion of the overdensities by the stronger radial 
drift, d) The maximum solid-to-gas ratio on the grid as a function of time. The average 
solids-to-gas ratio in the midplane is 0.5, whereas the maximum reaches well over ten 
times higher values in transient high pressure regions (yellow) and several hundred times 
higher values when the streaming instability is active (orange and blue). 
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Figure 2 Time series of the collapse of overdense seeds into gravitationally bound 
boulder clusters. The central panel shows the column densities of the four different sizes 
of boulders (in units of the mean column density of each size) plotted independently at 
a time just before self-gravity is turned on. All four particle sizes have concentrated at 
similar locations, an important prerequisite for the subsequent gravitational collapse. The 
surrounding panels show a time series of total column density of solids, in the radial- 
azimuthal {x-y) plane of the disc, summed over all particle sizes, starting from the upper 
left and progressing clockwise. Values are normalised to the average value across the 
grid (see colour bars in upper right panel). Times are given in orbital times T orb after self- 
gravity is turned on. Inset in each panel is an enlargement of a square region (indicated 
in the main panel) centred around the Hill sphere of the most massive cluster in the sim- 
ulation, represented by the white circle. These inserts show the log of the column density 
ratio (see colour bar in upper right panel) to capture the extreme values reached. Over- 
dense bands initially contract radially, forming thin filaments with densities high enough 
for a full non-axisymmetric collapse into gravitationally bound clumps to take place. As 
time progresses, the Hill sphere increases in radius as the clusters grow in mass by ac- 
creting boulders from the turbulent flow (see Supplementary Video for an animation of 
this simulation). 



15 




(this is Figure 2) 
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Figure 3 Mass accretion onto a gravitationally bound cluster. The plot shows the max- 
imum bulk density of solids as a function of time, normalised by the average gas density 
Drag force and vertical gravity are turned on at t = -10, while self-gravity and collisional 
cooling are turned on at t = 0. The density increases monotonically after the onset of 
self-gravity because gravitationally bound clusters of boulders form in the mid-plane. Af- 
ter only seven orbits peak densities in these clusters approach 10 4 p g or a million times 
the average boulder density in the disc. The coloured bars show the mass contained 
within the most massive Hill sphere in the box, in units of the mass of the 970 km ra- 
dius dwarf planet Ceres (M Ceres = 9.5 x l0 23 g). The most massive cluster accretes about 
0.5M Ce rcs per orbit (the entire box contains a total boulder mass of 50M Ce res)- The cluster 
consists of approximately equal fractions of the three larger boulder sizes. The smallest 
size, with O k t { = 0.25, is initially underrepresented with a fraction of only 15% because 
of the stronger aerodynamic coupling of those particles to the gas, but the fraction of 
small particles increases with time as the cluster grows massive enough to attract smaller 
particles as well. The mean free path inside the bound clusters is shorter than the size 
of the cluster, so any fragments formed in catastrophic collisions between the boulders 
will be swept up by the remaining boulders before being able to escape the cluster (see 
Supplementary Information). 
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